package ua.lviv.franko.integration;

public class TrapezoidalRuleQuadrature implements Quadrature {

	private int	numOfNodes;
	double		h;

	public TrapezoidalRuleQuadrature(int numberOfNodes) {
		this.numOfNodes = numberOfNodes;
		this.h = (2.0) / (double) numOfNodes;
	}

	@Override
	public double getW(int pos) {
		if ((pos == 0) || (pos == numOfNodes)) {
			return 0.5;
		} else {
			return 1.0;
		}
	}

	@Override
	public double getX(int pos) {

		return -1.0 + pos * h;
	}

	@Override
	public double integrate(double a, double b, IFunction1D func) {
		double res = 0.0;
		for (int i = 0; i <= numOfNodes; i++) {
			res += getW(i) * func.calculate((a + b) / 2 + (b - a) / 2 * getX(i));
		}
		res *= (b - a) / (double) numOfNodes;
		return res;
	}

	@Override
	public double[][] integrate(double a, double b, MatrixFunction1D func) {
		int matrixSize = func.matrixSize();

		double[][] res = new double[matrixSize][matrixSize];

		for (int i = 0; i <= numOfNodes; i++) {
			double[][] tmp = func.calculate((a + b) / 2 + (b - a) / 2 * getX(i));
			for (int j = 0; j < res.length; j++) {
				for (int k = 0; k < res.length; k++) {
					res[j][k] += getW(i) * tmp[j][k];
				}
			}
		}
		double aMinusB = (b - a) / (double) numOfNodes;
		for (int j = 0; j < res.length; j++) {
			for (int k = 0; k < res.length; k++) {
				res[j][k] *= aMinusB;
			}
		}

		return res;
	}

	@Override
	public double[][] integrate(double a, double b, MatrixFunction1DReqtangular func) {

		double[][] res = new double[func.matrixRows()][func.matrixColumns()];

		for (int i = 0; i <= numOfNodes; i++) {
			double[][] tmp = func.calculate((a + b) / 2 + (b - a) / 2 * getX(i));
			for (int j = 0; j < res.length; j++) {
				for (int k = 0; k < res[j].length; k++) {
					res[j][k] += getW(i) * tmp[j][k];
				}
			}
		}
		double aMinusB = (b - a) / (double) numOfNodes;
		for (int j = 0; j < res.length; j++) {
			for (int k = 0; k < res[j].length; k++) {
				res[j][k] *= aMinusB;
			}
		}

		return res;
	}

	@Override
	public double[] integrate(double a, double b, VectorFunction1D func) {
		int vecSize = func.vectorSize();

		double[] res = new double[vecSize];

		for (int i = 0; i <= numOfNodes; i++) {
			double[] tmp = func.calculate((a + b) / 2 + (b - a) / 2 * getX(i));
			for (int j = 0; j < res.length; j++) {
				res[j] += getW(i) * tmp[j];
			}
		}
		for (int j = 0; j < res.length; j++) {
			res[j] *= (b - a) / (double) numOfNodes;
		}

		return res;
	}

	public static void main(String[] args) {
		TrapezoidalRuleQuadrature g = new TrapezoidalRuleQuadrature(1000);
		System.err.println(g.integrate(0, 1, new IFunction1D() {

			@Override
			public double calculate(double x) {
				// TODO Auto-generated method stub
				return x * x * 3;
			}

			@Override
			public double calculateDerivaty(double x) {
				// TODO Auto-generated method stub
				return 0;
			}
		}));
	}

}
